require(estimatr)

#### data preparation ####
FtF.data <- read.csv("FtF_data.csv")
online.data <- read.csv("online_data.csv")

## young-respondent subset of the face-to-face data
FtF.young.data <- subset(FtF.data, age < 45)

## control-group subset
FtF.ctrl.data <- subset(FtF.data, treatment == 0)
FtF.young.ctrl.data <- subset(FtF.young.data, treatment == 0)
online.ctrl.data <- subset(online.data, treatment == 0)

# dichotomize social media use
FtF.data$social.media.d <- (FtF.data$social.media > 1) * 1
online.data$social.media.d <- (online.data$social.media > 1) * 1

# dichotomize popularity in neighborhood
FtF.data$neighborhood.d <- (FtF.data$neighborhood > 2) * 1
online.data$neighborhood.d <- (online.data$neighborhood > 2) * 1

## function for subgroup comparison of SDB
SDB.group.RC1.test <- function(y, x, data, treat, block, DQ) {
  temp.data.1 <- subset(data, data[, x] == 0)
  n.obs.1 <- nrow(temp.data.1)
  list.reg.1 <- lm_robust(formula(paste0(y, " ~ ", treat)), data = temp.data.1, 
                          fixed_effects = paste0("~ ", block))
  list.est.1 <- as.numeric(list.reg.1$coefficients)
  list.se.1 <- as.numeric(list.reg.1$std.error)
  DQ.est.1 <- mean(temp.data.1[temp.data.1[, treat] == 0, DQ], na.rm = TRUE)
  DQ.se.1 <- sqrt(DQ.est.1 * (1 - DQ.est.1) / 
                    sum(! is.na(temp.data.1[temp.data.1[, treat] == 0, DQ])))
  SDB.est.1 <- DQ.est.1 - list.est.1
  SDB.se.1 <- sqrt(list.se.1 ^ 2 + DQ.se.1 ^ 2)
  SDB.p.1 <- (1 - pnorm(abs(SDB.est.1 / SDB.se.1))) * 2
  temp.data.2 <- subset(data, data[, x] == 1)
  n.obs.2 <- nrow(temp.data.2)
  list.reg.2 <- lm_robust(formula(paste0(y, " ~ ", treat)), data = temp.data.2, 
                          fixed_effects = paste0("~ ", block))
  list.est.2 <- as.numeric(list.reg.2$coefficients)
  list.se.2 <- as.numeric(list.reg.2$std.error)
  DQ.est.2 <- mean(temp.data.2[temp.data.2[, treat] == 0, DQ], na.rm = TRUE)
  DQ.se.2 <- sqrt(DQ.est.2 * (1 - DQ.est.2) / 
                    sum(! is.na(temp.data.2[temp.data.2[, treat] == 0, DQ])))
  SDB.est.2 <- DQ.est.2 - list.est.2
  SDB.se.2 <- sqrt(list.se.2 ^ 2 + DQ.se.2 ^ 2)
  SDB.p.2 <- (1 - pnorm(abs(SDB.est.2 / SDB.se.2))) * 2
  SDB.diff.est <- SDB.est.2 - SDB.est.1
  SDB.diff.se <- sqrt(list.se.1 ^ 2 + DQ.se.1 ^ 2 + list.se.2 ^ 2 + DQ.se.2 ^ 2)
  SDB.diff.p <- (1 - pnorm(abs(SDB.diff.est / SDB.diff.se))) * 2
  out <- list(summary = c(SDB.est.1, SDB.est.2, SDB.diff.est, 
                          SDB.diff.est + qnorm(0.025) * SDB.diff.se, 
                          SDB.diff.est + qnorm(0.975) * SDB.diff.se, 
                          SDB.diff.p), 
              details = list(SDB.est.1 = SDB.est.1, SDB.se.1 = SDB.se.1, 
                             SDB.p.1 = SDB.p.1, n.obs.1 = n.obs.1, 
                             SDB.est.2 = SDB.est.2, SDB.se.2 = SDB.se.2, 
                             SDB.p.2 = SDB.p.2, n.obs.2 = n.obs.2, 
                             SDB.diff.est = SDB.diff.est, 
                             SDB.diff.se = SDB.diff.se, 
                             SDB.diff.p = SDB.diff.p))
  names(out$summary) <- c("SDB est (x = 0)", "SDB est (x = 1)", "SDB diff est", 
                          "SDB diff lower", "SDB diff upper", "SDB diff p")
  out
}

#### DQ-based estimates based only on the control group (Appendix B.1) ####
## test of Hypothesis 1 (FtF)
DQ.est.RC1.FtF <- mean(FtF.ctrl.data$DQ, na.rm = TRUE)
DQ.se.RC1.FtF <- sqrt(DQ.est.RC1.FtF * (1 - DQ.est.RC1.FtF) / 
                        sum(! is.na(FtF.ctrl.data$DQ)))

list.result.FtF <- lm_robust(list ~ treatment, data = FtF.data, 
                             fixed_effects = ~ block)
list.est.FtF <- as.numeric(list.result.FtF$coefficients)
list.se.FtF <- as.numeric(list.result.FtF$std.error)

SDB.est.RC1.FtF <- DQ.est.RC1.FtF - list.est.FtF
SDB.se.RC1.FtF <- sqrt(list.se.FtF ^ 2 + DQ.se.RC1.FtF ^ 2)

## test of Hypothesis 1 (online)
DQ.est.RC1.online <- mean(online.ctrl.data$DQ, na.rm = TRUE)
DQ.se.RC1.online <- sqrt(DQ.est.RC1.online * (1 - DQ.est.RC1.online) / 
                           sum(! is.na(online.ctrl.data$DQ)))

list.result.online <- lm_robust(list ~ treatment, data = online.data, 
                                fixed_effects = ~ block)
list.est.online <- as.numeric(list.result.online$coefficients)
list.se.online <- as.numeric(list.result.online$std.error)

SDB.est.RC1.online <- DQ.est.RC1.online - list.est.online
SDB.se.RC1.online <- sqrt(list.se.online ^ 2 + DQ.se.RC1.online ^ 2)

## test of Hypothesis 2
DQ.est.RC1.FtF.young <- mean(FtF.young.ctrl.data$DQ, na.rm = TRUE)
DQ.se.RC1.FtF.young <- sqrt(DQ.est.RC1.FtF.young * (1 - DQ.est.RC1.FtF.young) / 
                              sum(! is.na(FtF.young.ctrl.data$DQ)))

list.result.FtF.young <- lm_robust(list ~ treatment, data = FtF.young.data, 
                                   fixed_effects = ~ block)
list.est.FtF.young <- as.numeric(list.result.FtF.young$coefficients)
list.se.FtF.young <- as.numeric(list.result.FtF.young$std.error)

SDB.est.RC1.FtF.young <- DQ.est.RC1.FtF.young - list.est.FtF.young
SDB.se.RC1.FtF.young <- sqrt(list.se.FtF.young ^ 2 + DQ.se.RC1.FtF.young ^ 2)

SDB.diff.est.RC1 <- SDB.est.RC1.FtF.young - SDB.est.RC1.online
SDB.diff.se.RC1 <- sqrt(list.se.FtF.young ^ 2 + DQ.se.RC1.FtF.young ^ 2 + 
                          list.se.online ^ 2 + DQ.se.RC1.online ^ 2)

## tests of Hypotheses 3-6 (FtF)
EJK.RC1.result.FtF <- SDB.group.RC1.test("list", "EJK", FtF.data, 
                                         "treatment", "block", "DQ")

RT.RC1.result.FtF <- SDB.group.RC1.test("list", "RT", FtF.data, 
                                        "treatment", "block", "DQ")

social.media.RC1.result.FtF <- SDB.group.RC1.test("list", "social.media.d", 
                                                  FtF.data, 
                                                  "treatment", "block", "DQ")

neighborhood.RC1.result.FtF <- SDB.group.RC1.test("list", "neighborhood.d", 
                                                  FtF.data, 
                                                  "treatment", "block", "DQ")

## tests of Hypotheses 3-6 (online)
EJK.RC1.result.online <- SDB.group.RC1.test("list", "EJK", online.data, 
                                            "treatment", "block", "DQ")

RT.RC1.result.online <- SDB.group.RC1.test("list", "RT", online.data, 
                                           "treatment", "block", "DQ")

social.media.RC1.result.online <- SDB.group.RC1.test("list", "social.media.d", 
                                                     online.data, 
                                                     "treatment", "block", "DQ")

neighborhood.RC1.result.online <- SDB.group.RC1.test("list", "neighborhood.d", 
                                                     online.data, 
                                                     "treatment", "block", "DQ")

#### Figure A.2 ####
cairo_pdf("Figure_A2.pdf", width = 4.5, height = 3.2, pointsize = 9)
par(mar = c(4, 0, 1, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-0.9, 1), ylim = c(0, 14), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-0.05, c(1, 3:5, 7:9, 11:13), 1.05, c(1, 3:5, 7:9, 11:13), 
         lty = 3, col = "gray")
abline(v = seq(0, 1, 0.2), col = "gray")
points(DQ.est.RC1.FtF, 13, pch = 4)
segments(DQ.est.RC1.FtF + qnorm(0.025) * DQ.se.RC1.FtF, 13, 
         DQ.est.RC1.FtF + qnorm(0.975) * DQ.se.RC1.FtF, 13)
points(list.est.FtF, 12, pch = 4)
segments(list.est.FtF + qnorm(0.025) * list.se.FtF, 12, 
         list.est.FtF + qnorm(0.975) * list.se.FtF, 12)
points(SDB.est.RC1.FtF, 11, pch = 19)
segments(SDB.est.RC1.FtF + qnorm(0.025) * SDB.se.RC1.FtF, 11, 
         SDB.est.RC1.FtF + qnorm(0.975) * SDB.se.RC1.FtF, 11)
points(DQ.est.RC1.FtF.young, 9, pch = 4)
segments(DQ.est.RC1.FtF.young + qnorm(0.025) * DQ.se.RC1.FtF.young, 9, 
         DQ.est.RC1.FtF.young + qnorm(0.975) * DQ.se.RC1.FtF.young, 9)
points(list.est.FtF.young, 8, pch = 4)
segments(list.est.FtF.young + qnorm(0.025) * list.se.FtF.young, 8, 
         list.est.FtF.young + qnorm(0.975) * list.se.FtF.young, 8)
points(SDB.est.RC1.FtF.young, 7, pch = 19)
segments(SDB.est.RC1.FtF.young + qnorm(0.025) * SDB.se.RC1.FtF.young, 7, 
         SDB.est.RC1.FtF.young + qnorm(0.975) * SDB.se.RC1.FtF.young, 7)
points(DQ.est.RC1.online, 5, pch = 4)
segments(DQ.est.RC1.online + qnorm(0.025) * DQ.se.RC1.online, 5, 
         DQ.est.RC1.online + qnorm(0.975) * DQ.se.RC1.online, 5)
points(list.est.online, 4, pch = 4)
segments(list.est.online + qnorm(0.025) * list.se.online, 4, 
         list.est.online + qnorm(0.975) * list.se.online, 4)
points(SDB.est.RC1.online, 3, pch = 19)
segments(SDB.est.RC1.online + qnorm(0.025) * SDB.se.RC1.online, 3, 
         SDB.est.RC1.online + qnorm(0.975) * SDB.se.RC1.online, 3)
points(SDB.diff.est.RC1, 1, pch = 19)
segments(SDB.diff.est.RC1 + qnorm(0.025) * SDB.diff.se.RC1, 1, 
         SDB.diff.est.RC1 + qnorm(0.975) * SDB.diff.se.RC1, 1)
text(-0.9, 14, "Face-to-face survey (entire)", pos = 4, font = 2)
text(-0.8, 13, "Direct question", pos = 4)
text(-0.8, 12, "List question", pos = 4)
text(-0.8, 11, "Difference b/w questions", pos = 4)
text(-0.9, 10, "Face-to-face survey (18\u201344)", pos = 4, font = 2)
text(-0.8, 9, "Direct question", pos = 4)
text(-0.8, 8, "List question", pos = 4)
text(-0.8, 7, "Difference b/w questions", pos = 4)
text(-0.9, 6, "Online survey", pos = 4, font = 2)
text(-0.8, 5, "Direct question", pos = 4)
text(-0.8, 4, "List question", pos = 4)
text(-0.8, 3, "Difference b/w questions", pos = 4)
text(-0.9, 1, "Difference in differences\nb/w survey modes", pos = 4, font = 2)
axis(1, seq(0, 1, 0.2), lwd = 0.5)
mtext("Approval rate (or its difference)", side = 1, line = 2.5, at = 0.5)
dev.off()

#### Figure A.3 ####
cairo_pdf("Figure_A3.pdf", width = 5, height = 7, pointsize = 9)
layout(matrix(1:2, 2, 1))
par(mar = c(4, 0, 3, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-2.5, 1), ylim = c(0, 16), 
     main = "Face-to-face survey", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-0.8, c(1:3, 5:7, 9:11, 13:15), 1.05, c(1:3, 5:7, 9:11, 13:15), 
         lty = 3, col = "gray")
abline(v = seq(-0.75, 1, 0.25), col = "gray")
points(EJK.RC1.result.FtF$details$SDB.est.1, 15, pch = 4)
segments(EJK.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * EJK.RC1.result.FtF$details$SDB.se.1, 15, 
         EJK.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * EJK.RC1.result.FtF$details$SDB.se.1, 15)
points(EJK.RC1.result.FtF$details$SDB.est.2, 14, pch = 4)
segments(EJK.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * EJK.RC1.result.FtF$details$SDB.se.2, 14, 
         EJK.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * EJK.RC1.result.FtF$details$SDB.se.2, 14)
points(EJK.RC1.result.FtF$details$SDB.diff.est, 13, pch = 19)
segments(EJK.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * EJK.RC1.result.FtF$details$SDB.diff.se, 13, 
         EJK.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * EJK.RC1.result.FtF$details$SDB.diff.se, 13)
points(RT.RC1.result.FtF$details$SDB.est.1, 11, pch = 4)
segments(RT.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * RT.RC1.result.FtF$details$SDB.se.1, 11, 
         RT.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * RT.RC1.result.FtF$details$SDB.se.1, 11)
points(RT.RC1.result.FtF$details$SDB.est.2, 10, pch = 4)
segments(RT.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * RT.RC1.result.FtF$details$SDB.se.2, 10, 
         RT.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * RT.RC1.result.FtF$details$SDB.se.2, 10)
points(RT.RC1.result.FtF$details$SDB.diff.est, 9, pch = 19)
segments(RT.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * RT.RC1.result.FtF$details$SDB.diff.se, 9, 
         RT.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * RT.RC1.result.FtF$details$SDB.diff.se, 9)
points(social.media.RC1.result.FtF$details$SDB.est.1, 7, pch = 4)
segments(social.media.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * social.media.RC1.result.FtF$details$SDB.se.1, 7, 
         social.media.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * social.media.RC1.result.FtF$details$SDB.se.1, 7)
points(social.media.RC1.result.FtF$details$SDB.est.2, 6, pch = 4)
segments(social.media.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * social.media.RC1.result.FtF$details$SDB.se.2, 6, 
         social.media.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * social.media.RC1.result.FtF$details$SDB.se.2, 6)
points(social.media.RC1.result.FtF$details$SDB.diff.est, 5, pch = 19)
segments(social.media.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * social.media.RC1.result.FtF$details$SDB.diff.se, 5, 
         social.media.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * social.media.RC1.result.FtF$details$SDB.diff.se, 5)
points(neighborhood.RC1.result.FtF$details$SDB.est.1, 3, pch = 4)
segments(neighborhood.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * neighborhood.RC1.result.FtF$details$SDB.se.1, 3, 
         neighborhood.RC1.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * neighborhood.RC1.result.FtF$details$SDB.se.1, 3)
points(neighborhood.RC1.result.FtF$details$SDB.est.2, 2, pch = 4)
segments(neighborhood.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * neighborhood.RC1.result.FtF$details$SDB.se.2, 2, 
         neighborhood.RC1.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * neighborhood.RC1.result.FtF$details$SDB.se.2, 2)
points(neighborhood.RC1.result.FtF$details$SDB.diff.est, 1, pch = 19)
segments(neighborhood.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * neighborhood.RC1.result.FtF$details$SDB.diff.se, 1, 
         neighborhood.RC1.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * neighborhood.RC1.result.FtF$details$SDB.diff.se, 1)
text(-2.5, 16, "Extrajudicial killing", pos = 4, font = 2)
text(-2.35, 15, paste0("Do not know a victim (N=", 
                       format(EJK.RC1.result.FtF$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 14, paste0("Know a victim (N=", 
                       format(EJK.RC1.result.FtF$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 13, "Difference b/w groups", pos = 4)
text(-2.5, 12, "Red-tagging", pos = 4, font = 2)
text(-2.35, 11, paste0("Do not know a victim (N=", 
                       format(RT.RC1.result.FtF$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 10, paste0("Know a victim (N=", 
                       format(RT.RC1.result.FtF$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 9, "Difference b/w groups", pos = 4)
text(-2.5, 8, "Social media use", pos = 4, font = 2)
text(-2.35, 7, paste0("60 minutes or less (N=", 
                      format(social.media.RC1.result.FtF$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 6, paste0("More than 60 minutes (N=", 
                      format(social.media.RC1.result.FtF$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 5, "Difference b/w groups", pos = 4)
text(-2.5, 4, "Popularity in neighborhood", pos = 4, font = 2)
text(-2.35, 3, paste0("Dissatisfied (N=", 
                      format(neighborhood.RC1.result.FtF$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 2, paste0("Satisfied (N=", 
                      format(neighborhood.RC1.result.FtF$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 1, "Difference b/w groups", pos = 4)
axis(1, seq(-0.75, 1, 0.25), lwd = 0.5, labels = NA)
mtext(sprintf("%4.2f", seq(-0.75, 1, 0.25)), 
      side = 1, line = 1, at = seq(-0.75, 1, 0.25))
mtext("Magnitude of SDB (or its difference)", side = 1, line = 2.5, at = 0.125)
plot(NULL, NULL, bty = "n", xlim = c(-2.5, 1), ylim = c(0, 16), 
     main = "Online survey", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-0.8, c(1:3, 5:7, 9:11, 13:15), 1.05, c(1:3, 5:7, 9:11, 13:15), 
         lty = 3, col = "gray")
abline(v = seq(-0.75, 1, 0.25), col = "gray")
points(EJK.RC1.result.online$details$SDB.est.1, 15, pch = 4)
segments(EJK.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.025) * EJK.RC1.result.online$details$SDB.se.1, 15, 
         EJK.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.975) * EJK.RC1.result.online$details$SDB.se.1, 15)
points(EJK.RC1.result.online$details$SDB.est.2, 14, pch = 4)
segments(EJK.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.025) * EJK.RC1.result.online$details$SDB.se.2, 14, 
         EJK.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.975) * EJK.RC1.result.online$details$SDB.se.2, 14)
points(EJK.RC1.result.online$details$SDB.diff.est, 13, pch = 19)
segments(EJK.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.025) * EJK.RC1.result.online$details$SDB.diff.se, 13, 
         EJK.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.975) * EJK.RC1.result.online$details$SDB.diff.se, 13)
points(RT.RC1.result.online$details$SDB.est.1, 11, pch = 4)
segments(RT.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.025) * RT.RC1.result.online$details$SDB.se.1, 11, 
         RT.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.975) * RT.RC1.result.online$details$SDB.se.1, 11)
points(RT.RC1.result.online$details$SDB.est.2, 10, pch = 4)
segments(RT.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.025) * RT.RC1.result.online$details$SDB.se.2, 10, 
         RT.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.975) * RT.RC1.result.online$details$SDB.se.2, 10)
points(RT.RC1.result.online$details$SDB.diff.est, 9, pch = 19)
segments(RT.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.025) * RT.RC1.result.online$details$SDB.diff.se, 9, 
         RT.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.975) * RT.RC1.result.online$details$SDB.diff.se, 9)
points(social.media.RC1.result.online$details$SDB.est.1, 7, pch = 4)
segments(social.media.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.025) * social.media.RC1.result.online$details$SDB.se.1, 7, 
         social.media.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.975) * social.media.RC1.result.online$details$SDB.se.1, 7)
points(social.media.RC1.result.online$details$SDB.est.2, 6, pch = 4)
segments(social.media.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.025) * social.media.RC1.result.online$details$SDB.se.2, 6, 
         social.media.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.975) * social.media.RC1.result.online$details$SDB.se.2, 6)
points(social.media.RC1.result.online$details$SDB.diff.est, 5, pch = 19)
segments(social.media.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.025) * social.media.RC1.result.online$details$SDB.diff.se, 5, 
         social.media.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.975) * social.media.RC1.result.online$details$SDB.diff.se, 5)
points(neighborhood.RC1.result.online$details$SDB.est.1, 3, pch = 4)
segments(neighborhood.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.025) * neighborhood.RC1.result.online$details$SDB.se.1, 3, 
         neighborhood.RC1.result.online$details$SDB.est.1 + 
           qnorm(0.975) * neighborhood.RC1.result.online$details$SDB.se.1, 3)
points(neighborhood.RC1.result.online$details$SDB.est.2, 2, pch = 4)
segments(neighborhood.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.025) * neighborhood.RC1.result.online$details$SDB.se.2, 2, 
         neighborhood.RC1.result.online$details$SDB.est.2 + 
           qnorm(0.975) * neighborhood.RC1.result.online$details$SDB.se.2, 2)
points(neighborhood.RC1.result.online$details$SDB.diff.est, 1, pch = 19)
segments(neighborhood.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.025) * neighborhood.RC1.result.online$details$SDB.diff.se, 1, 
         neighborhood.RC1.result.online$details$SDB.diff.est + 
           qnorm(0.975) * neighborhood.RC1.result.online$details$SDB.diff.se, 1)
text(-2.5, 16, "Extrajudicial killing", pos = 4, font = 2)
text(-2.35, 15, paste0("Do not know a victim (N=", 
                       format(EJK.RC1.result.online$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 14, paste0("Know a victim (N=", 
                       format(EJK.RC1.result.online$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 13, "Difference b/w groups", pos = 4)
text(-2.5, 12, "Red-tagging", pos = 4, font = 2)
text(-2.35, 11, paste0("Do not know a victim (N=", 
                       format(RT.RC1.result.online$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 10, paste0("Know a victim (N=", 
                       format(RT.RC1.result.online$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 9, "Difference b/w groups", pos = 4)
text(-2.5, 8, "Social media use", pos = 4, font = 2)
text(-2.35, 7, paste0("60 minutes or less (N=", 
                      format(social.media.RC1.result.online$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 6, paste0("More than 60 minutes (N=", 
                      format(social.media.RC1.result.online$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 5, "Difference b/w groups", pos = 4)
text(-2.5, 4, "Popularity in neighborhood", pos = 4, font = 2)
text(-2.35, 3, paste0("Dissatisfied (N=", 
                      format(neighborhood.RC1.result.online$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 2, paste0("Satisfied (N=", 
                      format(neighborhood.RC1.result.online$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 1, "Difference b/w groups", pos = 4)
axis(1, seq(-0.75, 1, 0.25), lwd = 0.5, labels = NA)
mtext(sprintf("%4.2f", seq(-0.75, 1, 0.25)), 
      side = 1, line = 1, at = seq(-0.75, 1, 0.25))
mtext("Magnitude of SDB (or its difference)", side = 1, line = 2.5, at = 0.125)
dev.off()